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Phase transition of the Ising model is investigated on a planar lattice that has a fractal structure. 
On the lattice, the number of bonds that cross the border of a finite area is doubled when the linear 
size of the area is extended by a factor of four. The free energy and the spontaneous magnetization 
of the system are obtained by means of the higher-order tensor renormalization group method. The 
system exhibits the order-disorder phase transition, where the critical indices are different from that 
of the square-lattice Ising model. An exponential decay is observed in the density matrix spectrum 
even at the critical point. It is possible to interpret that the system is less entangled because of the 
fractal geometry. 
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I. INTRODUCTION 

Phase transition and critical phenomena have been one 
of the central issues in statistical analyses of condensed 
matter physicsi. When the second-order phase transi¬ 
tion is observed, thermodynamic functions, such as the 
free energy, the internal energy, and the magnetization, 
show non-trivial behavior around the transition temper¬ 
ature T— This critical singularity reflects the absence 
of any scale length at , and the power-law behavior 
of thermodynamic functions around the transition is ex¬ 
plained by the concept of the renormalization groupii^^— . 

Analytic investigation of the renormalization group 
flow in (/?^-model shows that the Ising model exhibits 
a phase transition when the lattice dimension is larger 
than one, which is the lower critical dimension^ii. In 
a certain sense, the one-dimensional Ising model shows 
rescaled critical phenomena around T = 0. When the 
lattice dimension is larger than four, which is the upper 
critical dimension, and provided that the system is uni¬ 
form, then the Ising model on regular lattices exhibits 
mean-field-like critical behavior. 

Compared with critical phenomena on regular lattices, 
much less is known on fractal lattices. Renormalization 
flow is investigated by Gefen et al.^“— where correspon¬ 
dence between lattice structure and the value of critical 
indices is not fully understood in a quantitative manner. 
For example, the Ising model on the Sierpinski gasket 
does not exhibit phase transition at any finite temper¬ 
ature, although the Hausdorff dimension of the lattice, 
djj = In3/In2 Ri 1.585, is larger than on o^^d^ . The 
absence of the phase transition could be explained by 
the fact that the number of interfaces, i.e. the outgoing 
bonds from a finite area, does not increase when the size 
of the area is doubled on the gasket. A non-trivial feature 
of this system is that there is a logarithmic scaling be¬ 
havior in the internal energy toward zero temperature!^. 
The effect of anisotropy has been considered recently!^. 
In case of the Ising model on the Sierpinski carpet, pres¬ 
ence of the phase transition is proved!^, and its critical 





FIG. 1. Composition of the fractal lattice. Upper left: a 
local vertex around an Ising spin shown by the empty dot. 
Middle: the basic cluster which contains Ni = 12 vertices. 
Lower right: the extended cluster which contains N 2 = 12^ 
vertices. In each step of the system extension, the linear size 
of the system increases by the factor of 4, where only 12 units 
are linked, and where 4 units at the corners are missing, if it 
is compared with a 4 by 4 square cluster. 

indices were roughly estimated by Monte Carlo simula- 
tionsii. It should be noted that it is not easy to collect 
sufficient number of data plots for finite-size scalingiS on 
such fractal lattices by means of Monte Carlo simula¬ 
tions, because of the exponential blow-up of the number 
of sites in a unit of fractal. 

In this article, we investigate the Ising model on a pla¬ 
nar fractal lattice, shown in Fig.[T] The lattice consists of 
vertices around the lattice points, which are denoted by 
the empty dots in the figure, where there are Ising spins. 
The whole lattice is constructed by recursive extension 
processes, where the linear size of the system increases 
by the factor of four in each step. If the lattice was a 
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regular square one, 4x4= 16 units are connected in the 
extension process, whereas only 12 units are connected on 
this fractal lattice; 4 units are missing in the corners. As 
a result, the number of sites contained in a cluster after 
n extensions is = 12", and the Hausdorff dimension 
of this lattice is = Inl2/ln4 si 1.792. The number 
of outgoing bonds from a cluster is only doubled in each 
extension process since the sites and the bonds at each 
corner are missing. If we evaluate the lattice dimension 
from the relation 

M = (1) 


1 /2 

where the r.h.s. takes the value (cosh2Ar) when a = 

a', and (cosh 2AT) when a ^ a', and where Eq. 0 

holds under the condition 

= (cosh27?)^^^ (5) 

The new parameter K is then expressed as follows 



Thus, if we introduce a factor 


between the linear dimension L and the number of out¬ 
going bonds M, we have d = 1.5, since M is proportional 
to \/L on the fractal. Remark that the value is different 
from djj ~ 1.792 

We report the critical behavior of the Ising model on 
the fractal lattice when the system size is large enough. 
Thermodynamic properties of the system are numerically 
studied by means of the Higher-Order Tensor Renormal¬ 
ization Group (HOTRG) methodic. The system exhibits 
the order-disorder phase transition, where the critical in¬ 
dices are different from the square lattice Ising model. 
In the next Section we introduce a representation of the 
Ising model in terms of a vertex model, which is suitable 
for numerical analyses by means of the HOTRG method. 
In Sec. HI, we show the calculated result around the tran¬ 
sition temperature . Gonclusions are summarized in 
the last Section. 


II. VERTEX REPRESENTATION 

We introduce a representation of the Ising model as 
a (symmetric) 16-vertex model. The Ising interaction 
between two adjacent Ising spins a and a', where each 
one takes either 4-1 or —1, is represented by the diagonal 
Hamiltonian 


H{cT^a') =—Jaa', (2) 

where J > 0 represents the ferromagnetic coupling. 
Throughout this article we assume that there is no exter¬ 
nal magnetic field. The corresponding local Boltzmann 
weight on the bond is given by 



i7(cr, cr') 


[ J /] 

exp 

kBT 

= exp 

-—— aa 

k^T 




(3) 


where T is the temperature, fee is the Boltzmann con¬ 
stant, and we have introduced a parameter K = J/k^T. 

It is possible to factorize the bond weight into 

two parts, by introducing an auxiliary spin s = ±1, which 
is often called an ‘ancilla’, and which is located between 
a and a—. A key relation is 


eKcrcr' 


1 

2 (cosh 2K) 




(4) 


= e 


K<ts 


2 (cosh 27?) 



(7) 


for each division of a bond, we can rewrite the Ising in¬ 
teraction in the following form 


( 8 ) 

S 

By means of the factorization in Eq. (|S]), we can map 
the square-lattice Ising model into the symmetric 16- 
vertex model, where the local vertex weight is defined 
as 


Tss's"s"' — ^ ^(TS 4Ecrs" IEcts"' • (9) 

<7 

In the upper-left corner of Fig. [U we have shown the 
graphical representation of the vertex weight 
where the open circle denotes the Ising spin a, which is 
summed up. The four short bars around the Ising spin 
in Fig. [T] show the halves of the bonds, where there are 
auxiliary spins s, s', s", and s'" at the end of each short 
bar. 

In case we consider a finite-size cluster with rectangular 
shape with free boundary conditions, we have to prepare 
a new boundary Boltzmann weight 


( 10 ) 

<7 

and a corner Boltzmann weight 


(7 


It should be noted that these boundary weights P^s's" 
and Cj, g, are obtained by taking partial trace for the 
vertex weight; we have the relations 


and 


y ,„T , „ n 
p _ ss s s 


y w 

L^s'" <7 
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( 12 ) 


(13) 
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where one can neglect the denominator when one is in¬ 
terested in the critical singularity; the denominators just 
subtract a regular function from the free energy of the 
system. In case that one needs fixed boundary condi¬ 
tions, it is sufficient to avoid taking the configuration 
sum for (T in the r.h.s. of both Eq. m and Eq. (HU, 
and to set all the boundary spins to be either -|-1 or — 1 
according to the condition. The vertex weights 
-f’ss's"! ^ss' invariant under arbitrary permuta¬ 
tion of the indices. 

There are various choices of the factorization of the 
bond weight in Eq. (|S]). Instead of using the relation in 
Eq. ([7]) , one can introduce an asymmetric decomposition 


W = 


V cosh K, 
\/coshK, 


Vsinh K \ 

—\/sinhitr J ’ 


(14) 


where we have used the matrix notation for the weight 
This expression is often employed in the tensor 
network formulations^^, which does not require any typ¬ 
ical symmetry for local weights, as long as the numerical 
treatment is concerned. In case this asymmetric factor¬ 
ization in Eq. (ITTl) is employed, one has to care about the 
order of the indices in In the following nu¬ 

merical calculation, we use the symmetric factorization. 

The fractal lattice we treat in this article is constructed 
by a recursive joining process of the local tensors, which 
is nothing but a vertex weight in Eq. dH]) at the begin¬ 
ning. In each extension process, we join 12 local tensors 
as shown in the middle of Fig. [TJ In the joining process, 
we take the configuration sum for those tensor indices in¬ 
side the cluster, leaving those on the border that become 
new tensor indices of the extended tensor. Because of 
the fractal geometry, some of the bonds inside the cluster 
are not connected with each other. We also take config¬ 
uration sum for these dangling bonds, and the process 
just change the normalization of the partition function 
by amount of 


boundary conditions to evaluate 

ZJT) = , (16) 


where is the renormalized local tensor obtained after 
n extensions. 


III. NUMERICAL RESULTS 

In order to simplify the numerical analysis, we choose 
the parameterization J = A:b = 1, and thus we have 
K = 1/T. In the numerical calculation by means of 
HOTRG, we keep I? = 24 states at most for block spin 
variables. We have verified that the choice £) = 24 is 
sufficient for obtaining the converged free energy 

F^{T) = -kBT In ZJT) (17) 

in the entire temperature region^^. We treat the free 
energy per site 


F (T) 

f{T) = hm ^ (18) 

n-^-oo 

in the following thermodynamic analyses, where the r.h.s. 
converges already for n 30 . 

Figure [5] shows the temperature dependence of the spe¬ 
cific heat per site 


c{T) = ^u{T), (19) 

where u(T) is the internal energy per site 

<T) = A , (20) 




2 cosh K 


2 (cosh 2 A) 



(15) 


for each, if we choose the definition of in Eq. 0. 
We take the rescaling effect into account, although the 
rescaling is not essential to the thermodynamic properties 
of the system, in particular to its critical singularity. In 
this manner, what we are dealing with is the Ising model, 
where there are only spins denoted by the empty dots in 
Fig.dl 

At first we have only 4 spins s, s', s", and s'" on the 
outgoing bonds, and after n extensions of the system, we 
have 4x2" border spins on the surface of the extended 
cluster. The application of the HOTRG to this fractal 
system is straightforward. The recursive structure of the 
lattice is suitable for the repeated process of system ex¬ 
tensions and renormalization group transformations in 
the HOTRG method. The partition function ZJT) of 
the system after n extensions is obtained by the con¬ 
traction of the extended tensors; we choose the periodic 


and the temperature derivatives are performed numeri¬ 
cally. There is no singularity in c{T) around its maxi¬ 
mum. One might find a weak non-analytic behavior at 
Tj. Ri 1.317, which is marked by the dotted line in the 
figure; the numerical derivative of c(T) with respect to 
temperature (plotted in the inset) has a sharp peak at 
the critical temperature T^. It is, however, difficult to 
determine the critical exponent a precisely, because of 
the weakness in the singularity; as shown in the figure, 
c(T) around T^ is almost linear in T, and therefore a is 
nearly zero. 

Figure [3] shows the spontaneous magnetization per site 
m{T), which is obtained by inserting a cr-dependent local 
weight 


Ts s's"s'" = E (21) 

(7 

into the system. Since the fractal lattice is inhomoge¬ 
neous, the value is weakly dependent on the location of 
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FIG. 2. The specific heat c{T) per site in Eq. (1191) . Inset: 
numerical derivative of c{T) with respect to temperature; a 
sharp peak is observed at « 1.317. 



T 


FIG. 3. The spontaneous magnetization per site m{T). Inset: 
the power-law behavior below = 1.31716. 

the observation site, but the critical behavior is not af¬ 
fected by the location; we choose a site from the four 
sites that are in the middle of the 12-site cluster shown 
in Fig.[T] The numerical calculation by HOTRG captures 
the spontaneous magnetization m(T) below since any 
tiny round-off error is sufficient for breaking the sym¬ 
metry inside low-temperature ordered state. Around 
the transition temperature, the magnetization satisfies 
a power-law behavior 

m(T) oc |T, (22) 

where the precision of the exponent is around 2%, which 
can be read out from the inset of Fig.Oas a tiny deviation 
from the linear dependence (the dashed lines) in to(T)^/^ 
near . 

As a byproduct of the numerical HOTRG calculation, 
we can roughly observe the entanglement spectrum,™ 
which is the distribution of the eigenvalue of the den¬ 
sity matrix that is created for the purpose of obtaining 



FIG. 4. Decay of the singular values after n = 8 extensions. 


the block spin transformation. Since the effect of envi¬ 
ronment is not considered in our implementation of the 
HOTRG method, the eigenvalue uj^ = Af is obtained as 
the square of the singular values Aj in the higher-order 
singular value decomposition applied to the extended ten¬ 
sors. Figure m shows Wj at T = in the decreasing or¬ 
der. The decay is rapid, and therefore further increase of 
the number of block-spin state from H = 24 to a larger 
number does not significantly improve the precision in 
the difference in f{T ^) between D = 8 and H = 16 is 
already of the order of 10“®. It should be noted that the 
eigenvalues are not distributed equidistantly in logarith¬ 
mic scale; the corner double line structure is absent^iiS^. 


IV. CONCLUSIONS AND DISCUSSIONS 

We have investigated the Ising model on the fractal 
lattice shown in Fig.[I]by means of the HOTRG method. 
The calculated specific heat c(T) suggests that the model 
shows 2'^'^ order phase transition. Qualitatively speak¬ 
ing, the presence of weak singularity in the specific heat 
agrees with the result of the e-expansion, which shows 
the increasing nature of the critical exponent in c{T) with 
respect to the space dimension d The calculated spon¬ 
taneous magnetization m{T) also supports the 2'^'^ or¬ 
der phase transition with the exponent /^fractal ~ 0.0137, 
which is smaller by one order of magnitude than the crit¬ 
ical exponent /Sgquare = 1/8 = 0.125 of the square-lattice 
Ising model. 

The fractal structure of the lattice modifies the en¬ 
tanglement spectrum from that on the square lattice ex¬ 
plained by the corner double line picture^iiS^. Since each 
corner is missing in the fractal structure in Fig. [U short- 
range entanglement is almost filtered out in the process of 
the renormalization group transformation. This may be 
the reason why we do not need many degrees of freedom 
for the renormalized tensors. The situation is similar to 
the entanglement structure reported in the tensor net- 
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work renormalization^^— . 

The lattice geometry of the fractal lattice can be mod¬ 
ified in several manners. For example, one can alternate 
the system extension process of the fractal for the pur¬ 
pose of modifying the Hausdorff dimension; for every odd 
n the extension with 12 vertices shown in Fig. [T] is per¬ 
formed, and for even n the normal extension with 16 ver¬ 
tices on the square-lattice is performed. Alternatively, 
one can also modify the basic cluster, in such a manner 
as introducing 6 by 6 cluster where 4 corners are missing, 
etc. It is also worth considering three-dimensional fractal 
lattice, and apply the HOTRG method as it was done for 
the cubic lattice Ising model^. These modifications do 
not spoil the applicability of the HOTRG method while 


the numerical requirement is heavier than the current 
research. Analyses of quantum systems on a variety of 
fractal lattice is another possible extensions^^^. These 
further study may clarify the role of the entanglement in 
the universality of the phase transition in both regular 
and fractal lattices. 
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